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ABSTRACT 

We use a population synthesis approach to characterise, as a function of cosmic time, 
the extragalactic close binary population descended from stars of low to intermediate 
initial mass. The unresolved gravitational wave (GW) background due to these systems 
is calculated for the 0.1-10 mHz frequency band of the planned Laser Interferometer 
Space Antenna (LISA). This background is found to be dominated by emission from 
close white dwarf-white dwarf pairs. The spectral shape can be understood in terms 
of some simple analytic arguments. To quantify the astrophysical uncertainties, we 
construct a range of evolutionary models which produce populations consistent with 
Galactic observations of close WD-WD binaries. The models differ in binary evolu- 
tion prescriptions as well as initial parameter distributions and cosmic star formation 
histories. We compare the resulting background spectra, whose shapes are found to 
be insensitive to the model chosen, and different to those found recently by Schneider 
et al. (2001). From this set of models, we constrain the amplitude of the extragalactic 
background to be 1 x 10~ 12 < figw(l mHz) < 6 x 10 -12 , in terms of f2 gw (/), the fraction 
of closure density received in gravitational waves in the logarithmic frequency interval 
around /. 
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1 INTRODUCTION 

Except at very low radio frequencies, most electromag- 
netic telescopes have good angular rejection, so that faint 
sources and backgrounds can be seen by looking between 
bright sources. In contrast all currently implemented grav- 
itational wave detectors, and most of those envisaged for 
the future, simultaneously respond to sources all over 
the sky, modified only by a beam pattern of typically 
quadrupole form. It is therefore important to understand 
the brightness of the gravitational wave sky, since this 
will limit the ultimate sensitivity attainable in gravita- 
tional wave astronomy. One immediate pressure to under- 
stand this background comes from the need to set de- 
sign requirements for the ESA/NASA Laser Interferometer 
Space Antenna (LISA) mission (LISA mission documents 
and status may be found at http://lisa.jpl.nasa.gov/ and 
http:/ /sci.esa.int/home/lisa/). 

In this paper, we attempt to predict the gravitational 
wave background produced by all the binary stars in the uni- 
verse, excluding neutron stars and black holes. This is be- 
lieved to be the principal source of gravitational wave back- 
ground in the frequency range 1(T 5 < / < 10~ 1 Hz. Below 
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10" 5 Hz, the background is probably dominated by merging 
supermassive black holes, and above 10" 1 Hz, it is proba- 
bly dominated by merging neutron stars and stellar mass 
black holes (whose complicated and poorly understood for- 
mation histories and birth velocities make predictions more 
uncertain, cf. Belczynski, Kalogera & Bulik 2002). 

Besides the extragalactic background, there is also a 
Galactic background produced by the binary stars in our 
Milky Way (Evans, Iben & Smarr 1987; Hils, Bender & Web- 
bink 1990; Nelemans, Yungelson & Portegies Zwart 2001c). 
Although the Galactic background is many times larger 
in amplitude than both the extragalactic background and 
LISA's design sensitivity, the individual binaries contribut- 
ing to it can be (spectrally) resolved and removed at fre- 
quencies above ~3x 10 -3 Hz (Cornish & Larson 2002). 
Below this frequency they cannot be removed (at least in 
a mission of reasonable lifetime ~ 3 years), but the unre- 
solved Galactic background will be quite anisotropic. As the 
detector beam pattern rotates about the sky, the Galactic 
background will thus be modulated, while the isotropic (or 
nearly so; see Kosenko & Postnov 2000) distant extragalac- 
tic background will not. Modelling of the angular distribu- 
tion of the Galactic background using both a priori models 
and the observed distribution of higher frequency resolved 
sources will thus allow the Galactic background to be sub- 
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tracted to some precision (Giampieri & Polnarev 1997). In 
addition, there will be anisotropics due to the distribution 
of local galaxies, at the level of 10 per cent of the distant ex- 
tragalactic background from the LMC, and at the per cent 
level from M31 or the Virgo cluster (see also Lipunov et al. 
1995). 

The immediate motivation for this work is a de- 
sign issue for LISA. One of LISA's major science 
goals (see the LISA Science Requirements document at 
http://www.tapir.caltech.edu/listwgl/) is the detection of 
gravitational waves from compact objects spiralling into su- 
permassive black holes (Finn & Thorne 2000; Hils & Bender 
1995), since these can provide precision tests of strong field 
relativity and the no-hair theorem (Hughes 2001). However, 
these signals are weak, and their templates not yet fully un- 
derstood. It has thus been proposed that LISA should be 
designed with somewhat greater sensitivity to increase the 
probability that these signals are detected. However, this 
would be pointless if the principal background were cosmo- 
logical rather than instrumental. As we shall see (Fig. 16), 
we find that this is most probably almost, but not quite the 
case at the relevant frequencies (4—10 mHz). So there would 
be a point to increasing LISA's sensitivity in the 4 — 10 mHz 
range, but not to increasing it by more than a factor of 3 in 
gravitational wave amplitude h (9 in Q oc f 2 h 2 ). 

A second motivation for this work comes from the 
fact that this background is an astrophysical foreground to 
searches (both with LISA and with future detectors with 
extended frequency range and sensitivity) for backgrounds 
produced in the very early universe. Gravitational waves 
from bubble walls and turbulence following the electroweak 
phase transition are expected to be in the LISA frequency 
band, with amplitude that could be well above LISA instru- 
mental sensitivity (Kamionkowski, Kosowski & Turner 1994; 
Kosowsky, Mack & Kahniashvili 2002; Apreda et al. 2002). 
Another potential source of isotropic gravitational waves in 
the LISA band are those produced when dimensions beyond 
the familiar four compactified, which occurred when the uni- 
verse had temperature kT >TeV (Hogan 2000). 

Note that detection of a gravitational wave background 
can possibly be made even if it is considerably below the 
noise limit of the LISA detectors shown in our Fig. 16. 
This can be done by comparing the signals from Michel- 
son beam combinations (sensitive to instrument noise and 
gravitational waves) with Sagnac beam combinations (sen- 
sitive to instrument noise, but insensitive to gravitational 
waves), thus calibrating the instrumental noise — cf. Tinto, 
Armstrong & Estabrook (2001), Hogan & Bender (2001). 

Gravitational waves are the only directly detectable 
relic of inflation in the early universe, and their detec- 
tion over a range of frequencies would provide a valuable 
test of models of inflation (Turner 1997). It has been pro- 
posed that advanced space-based gravitational wave detec- 
tors might search for the background of gravitational waves 
from inflation. The gravitational waves from slow-roll infla- 
tion models contribute to the critical density in the universe 
l!gw < 10~ 15 per octave of frequency. We shall see (Fig. 
8, 17) that the gravitational wave background from cosmo- 
logical binaries makes such detection impractical except at 
frequencies below 10 -5 Hz (where supermassive black holes 
continue to make it impossible), or above 0.1 Hz. 

A third motivation is that a detection of the extragalac- 



tic binary background, e.g. by LISA, would set an indepen- 
dent (and unaffected by dust extinction) constraint on a 
combination of the star formation history of the universe 
and binary star evolution. 

There have been previous estimates of the extragalactic 
binary background. Hils, Bender & Webbink (1990) made 
detailed estimates of the Galactic binary background, and 
estimated that the extragalactic background from close dou- 
ble white dwarf pairs should be about 2 per cent (in flux or Q 
units) of the Galactic background. This estimate was refined, 
using more modern star formation histories, by Kosenko & 
Postnov (1998), who found instead a level of ~10 per cent. 
Schneider et al. (2001) used a descendant of the Utrecht pop- 
ulation synthesis code to estimate the extragalactic binary 
background as a function of frequency, and claimed that the 
background should have a large peak at ~ 3 x 10~ 5 Hz, just 
below the frequency at which typical binaries have a lifetime 
that equals the age of the Universe. 

We have followed the spirit of this previous work, but 
with an independent binary population synthesis code. More 
importantly, we have devoted much effort to the normalisa- 
tion of the background, to understanding the contributions 
of different types of binaries and their formation pathways 
to the background, and to estimating the uncertainties in 
all of these, so that we can have a better idea of the sources 
and level of uncertainty in the predicted background. 

The paper is organised as follows: In section 2, we de- 
scribe the gravitational wave (GW) emission from a binary 
system, then in section 3 we outline the main evolutionary 
pathways to the close double degenerate (DD) stage, which 
we shall see is the dominant source of GW background in 
the LISA band. In section 4, we use the preceding sections to 
make some simple analytic arguments about the nature of 
DD inspiral spectra. We describe the use of the bse code 
in our population synthesis, in section 5, then go on to 
construct a set of synthesis models whose results we test 
against the observed Galactic DD population. We also mo- 
tivate some modifications made to the prescription for the 
evolution of AM CVn stars in the BSE code. In section 6, we 
present the cosmological integrals used in the code, along 
with the cosmic star formation history and overall normali- 
sation chosen. Section 7 is devoted to a discussion of the GW 
background spectra produced by our code, in terms of the 
systems contributing to the background and the progenitors 
of these sources. We also discuss the differences between our 
population synthesis models. In section 8, we place limits 
on the maximum and minimum expected background sig- 
nals, and compare these with the LISA sensitivity and in 
section 9 with previous work. In section 10 we summarise 
and conclude. 



2 GRAVITATIONAL WAVES FROM A 
BINARY SYSTEM 

A binary system of stars in circular orbit with masses Mi 
and M2 and orbital separation a emits gravitational radia- 
tion, at the expense of its orbital energy, at a rate given by 
(Peters & Mathews 1963) 

32 G 4 (MiM 2 ) 2 (Mi +M 2 ) 
5 c 5 a° 
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l.O x 10" 



.(M{M^ 2 (Mj+Mi) 
(a') 5 



erg s 



(1) 



where primes denote quantities expressed in solar units, i.e. 
M/M©, a/R©. The gravitational radiation is emitted at 
twice the orbital frequency v of the binary, / c i rc = 2v — fl/n. 

If the binary is eccentric with eccentricity e, this ex- 
pression must be generalised to include emission at all har- 
monics n of the orbital frequency, f n =w — nQ/2n, where 
A = (a _3 G(Mi +M 2 )) 1/2 . The luminosity in each harmonic 
is given by 



L(n,e) =g(n,e)L c 



(2) 



where L c irc is the luminosity of a circular binary with sep- 
aration a, as given in Eq. 1, where a is now the relative 
semi-major axis of the eccentric orbit, and the g(n, e) are 
defined in eq. (20) of Peters & Mathews (1963). The total 
specific luminosity Lf = dL(f)/df of the system is then a 
sum over all harmonics: 



L f( e ) = ^circ ^2 g{n,e)5(f - nv). 



The total luminosity is 

OQ 

L = Lcirc ^2 9( n , e ) = 



1 + 23 e 2 + 32 e 4 
'24, ~ 96 
(l_ e 2)7/2 ' 



(3) 



(4) 



For eccentric orbits, the emission spectrum of Eq. 3, 
g(n, e)L c i rc as a function of / = nv consists of points along 
a skewed bell-shaped curve with maximum near the relative 
angular velocity at pericentre, where the greatest acceler- 
ations are experienced (2nv ~ Q p , where Q p is the angu- 
lar velocity of the relative orbit at pericentre, v p /r p ). In 
terms of harmonic number, a good approximation for all 
e (becoming very good for e > 0.5) is that Lf peaks at 
n = 1.63(1 - e)~ 3/2 , and fL f peaks at n = 2.16(1 - e)- 3/2 . 



3 EVOLUTION TO THE DD STAGE 

We shall see that the GW background is dominated by the 
emission from close double degenerate (DD) binaries at fre- 
quencies 10~ 4 < / < 10 _1 Hz. In this work, the term DD will 
refer to WD-WD pairs and loosely to WD-naked helium 
star pairs, i.e. we exclude neutron stars from our definition. 
In this section we describe the two main evolutionary path- 
ways from the zero-age main sequence (ZAMS) to the close 
DD stage. The route followed depends mainly on the initial 
orbital separation of the ZAMS stars. Similar descriptions 
can be found in e.g. Webbink & Han (1998). 

We begin with an intermediate-mass ZAMS binary sys- 
tem with primary mass Mi, secondary mass M2 (< Mi), 
semi-major axis a and eccentricity e. The orbit may evolve 
somewhat due to tidal interactions between the stars, partic- 
ularly if they have convective envelopes. When the primary 
evolves off the main sequence and swells in size, it may fill 
its Roche lobe and start to transfer matter on to the sec- 
ondary. The stability of this mass transfer determines which 
of the two main pathways to the DD stage is commenced. 



3.1 CEE+CEE 

If the primary fills its Roche lobe when it has a deep 
convective envelope (i.e. on the red giant branch (RGB) 



or asymptotic giant branch (AGB)), then for mass ratios 
M1/M2 >0.6, the ensuing mass transfer is dynamically un- 
stable (for conservative transfer). The envelope of the pri- 
mary spills on to the secondary on a dynamical timescale, 
leading to the formation of a common envelope, inside which 
orbit the secondary and the core of the primary. The enve- 
lope is frictionally heated at the expense of the stars' orbital 
energy, until eventually either they coalesce, or the envelope 
is heated sufficiently that it is ejected from the system, leav- 
ing the primary's core (a hot subdwarf which will rapidly 
cool to become a WD, or if the primary was on the RGB 
and had mass Mi >2 M©, then a helium star which will 
evolve to the WD stage). The basic idea of the common 
envelope phase is well accepted and observationally moti- 
vated, though not well simulated (see e.g. Livio & Soker 
1988; Iben & Livio 1993; Taam & Sandquist 2000). Several 
formalisms have been proposed to model it in population 
synthesis studies. The evolution code used here (see section 
5.1) follows closely the prescription of Tout et al. (1997) 
(originally from Webbink 1984), in which 



-E?bind,i = a(E OI b,l — -Eorb.i), 



(5) 



where i?bind,i is the initial binding energy of the envelope 
of the overflowing giant star (or the sum of both envelopes' 
binding energies if both stars are giants), parametrized by 
-Ebind.i = — G/A(MiM en v,i/-Ri), where A is of order unity, 
and is calculated in the bse code (see section 5.3). -E or b,i 
and -E rb,f are respectively the initial and final orbital bind- 
ing energies of the core-plus-secondary system, and a is the 
so-called common envelope efficiency parameter, also of or- 
der unity, usually taken to be a parameter to be fitted to 
observations. Variations to this prescription will be consid- 
ered in sections 5.2 and 5.3. 

Continuing with the system's evolution, the secondary 
star later evolves off the main sequence, and a second com- 
mon envelope phase is likely to occur, leading to further 
orbital shrinkage. If once again the stars do not coalesce 
then we will be left with a close(r) pair of remnants, one 
or both of which may be helium stars, which in time will 
evolve to the WD stage. (It is not uncommon for either he- 
lium star to overflow its Roche lobe upon leaving the helium 
main sequence; this can lead to either stable mass transfer 
or to a futher common envelope phase.) In this picture, the 
second-formed WD will be the less massive of the pair, since 
the giant star from which it descended had a smaller core 
mass when its core growth was halted as it lost its envelope. 

3.2 Stable RLOF+CEE 

If Roche lobe overflow occurs when the primary is in the 
Hertzsprung gap, that is after the primary has exhausted 
its core hydrogen and before it has developed a deep con- 
vective envelope and ascended the giant branch, then Roche 
lobe overflow may be dynamically stable for moderate mass 
ratios, and a phase of stable but rapid mass transfer can 
occur. In this way, the primary transfers its envelope to the 
secondary, leaving a compact remnant, and a common enve- 
lope phase is avoided, since by the time the primary evolves 
to the giant branch, the mass ratio has been sufficiently in- 
verted that mass transfer remains dynamically stable. The 
orbital separation will typically have increased during this 
phase (for conservative mass transfer at least), since much of 
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the transfer was from the less-massive to the more-massive 
star. When the secondary evolves off the main sequence, it 
will most likely fill its Roche lobe on the RGB, so that a 
common envelope phase ensues, and a close DD is born, 
provided that the resulting orbital shrinkage does not lead 
to coalescence. The second-formed WD will this time be the 
more massive, since its progenitor was the more evolved at 
the time of its overflow. 

The initial conditions for this route occupy a smaller 
range in initial orbital semimajor axis than the CEE+CEE 
route, but as it results in the injection of DD systems only 
at very short periods, we expect both pathways to be sig- 
nificant contributors to the close DD population, i.e. those 
systems contributing to the GW background in the LISA 
waveband. We note also that both routes ought to lead to 
the production of DDs with circular orbits, even if the ZAMS 
eccentricity was non-zero, since tidal circularisation is rapid 
when a system contains a near-Roche lobe-filling convective 
star. 



4 ANALYTIC ARGUMENTS ABOUT 
SPECTRAL SHAPE 

Given only the above, we can make some predictions as to 
the shape of the GW spectrum seen today. A somewhat 
analogous treatment is given in Hils et al. (1990). We con- 
sider the evolution under GW emission of a population of 
DDs after creation as in Section 3, with circular orbits. We 
deal here with detached systems; the spectral shape due to 
interacting pairs is discussed in section 7.1.2. 

Here and throughout, we use v for orbital frequencies 
and / for gravitational wave frequencies. For circular orbits, 
f = 2v. 

The number density N(v,t) of binary WDs per unit 
orbital frequency interval at time t must obey the continuity 
equation 



for Eq. 6) Na{v,t; v 1 ,t') as a function of time is given by 



(6) 



where N b (v,t) is the birth rate (after nuclear evolution and 
mass transfer) of WD-WD systems per unit frequency. Now 
for a given source, we know that E OI b = — L gw , and using 
Eq. I along with E olh = -GM 1 M 2 /(2a) and Kepler's law, 
we obtain 



/3 



96 8/3 (GM\ 5/3 n 
= Kv 11/3 

~ (3.7 x I0-V 2 )(X/M o ) 5/ V 1/3 , (7) 
where we have used the definition of the chirp mass M, 



M 



M 3/5 M 3/5 

(Mi + M2) 1 / 5 



(8) 



We solve Eq. 7 to give the evolution v(t) for iVb = 
S(t — t',v — v'), i.e. for a single source injected at frequency 
v' at time t' 



v(t) 



-8/3 



.'-8/3 



8K(t' - t)/3. 



(9) 



Ng dv oc dt(v) 



dv 



_A45/3„ll/3 



8 t 



t' + 



8K 



(10) 



The corresponding source number density (Green's function 



since, as the system traces out a path in v, it spends a time 
at each point inversely proportional to its velocity v through 
frequency space. 

We then consider a real injection spectrum N b (v',t'), 
for lAnm < v' < lAnax- The resulting number density N(v,t) 
is given by 

N(u,t) = f f N b (v',t')N G (v,t;v',t')dt' dv'. (11) 
ii/ min Jo 

Since L gw oc v 10/3 M 10/3 , we can then construct the GW 
emission spectrum by taking F gw (/,t) oc f 10/3 M 10/3 N(f,t). 

The choice of DD injection spectrum is therefore in- 
strumental in determining the shape of the GW emission 
spectrum. We can estimate its shape as follows: we will later 
choose to distribute ZAMS orbital semimajor axis uniformly 
in log a, i.e. also uniformly in log v, for given initial M\ and 
M2. We suppose that, for at least the CEE+CEE route (see 
section 3) , the common envelope phases lead to some mean 
orbital shrinkage factor, so that WD-WD pairs at their birth 
are also distributed roughly uniformly in log v. We then have 
Nb(v') oc 1/v', from some z/ m i n <C v of interest, up to i/ mas 
(see also fig. 1 of Webbink & Han 1998). This is the max- 
imum orbital frequency at which a system can exit a com- 
mon envelope phase and survive to become a WD-WD pair. 
Upon CE exit, the newly exposed stellar core will be a hot 
subdwarf, larger than the WD it will cool to become, or it 
could be a naked helium star, which will eventually evolve to 
the WD stage. The maximum injection frequency at WD- 
WD birth is set by the minimum orbital separation that will 
keep this object (and the first-formed WD) from overflowing 
its Roche lobe on the way to the WD stage, whether this is 
at the exit of common envelope or (applicable to the helium 
star case) as its radius changes due to nuclear evolution. 

For illustration, we compute the emergent spectrum for 
a fiducial population of 0.5 M Q WD-WD pairs. The radius 
of a 0.5 Mq naked helium star does not exceed ~ 0.13 Rq 
on its way to the WD stage, which sets i/ m « ~ 0.7 mHz. 
If we then assign a constant pair formation rate, so that 
Nb(v',t') = Nb{v'), and perform the integral in Eq. 11, wc 
obtain the spectral shape shown in Fig. 1. Note that the 
spectrum is truncated at a frequency above which the inspi- 
ralling WDs would undergo Roche lobe overflow and merge, 

/merge — 2l^ m crgc — 40 HlHz. 

If instead we only inject sources for < t < r, and look 
at the spectrum obtained for t > r = 1 Gyr, (Fig. 2) , we see 
that the basic spectral shape is little affected. 

Because of the strong dependence of v on v, a given 
system of specified age will either have merged or will have 
remained at essentially constant separation. Thus there are 
two clear physical regimes displayed in the spectra, sepa- 
rated by the injection frequency from which a source could 
have reached contact due to GW losses in the time t since 
its birth, v ciit ~ 0.03 - 0.04 mHz for t ~ 5 - 10 Gyr. (In all 
relevant situations for us, v CI - lt < fmax-) 

At / < 2^ cr it lies a 'static regime', in which losses due to 
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Figure 1. Gravitational wave spectrum arising from constant 
WD-WD formation rate, at times 2, 5, 10, 100 and 1000 Gyr, 
increasing in the direction of the arrow shown. 
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Figure 2. Gravitational wave spectrum arising from a burst of 
WD— WD formation between and 1 Gyr. Curves plotted are 
spectra at times 2, 5, 10, 100 and 1000 Gyr, increasing in the 
direction of the arrow shown. 



GW are negligible in the time available, giving N{v) oc v 1 
and hence /F gw oc f w f 3 M 10/z . For />2z/ cri t, we are in the 
'spiral-in' regime. In the case of a burst of DD formation 
(Fig. 2), sources simply sweep through this region on the 
way to merger, so that we have N(u) oc v^ 11 ^ 3 M~ 5 / 3 , giving 
/F gw oc f 2 / 3 M 5/3 . If we have a constant DD formation 
rate (Fig. 1), then for 2^ cr it < / < 2f max , merging systems 
are continually being injected, so that N(u) is less steeply 
decreasing than v~ lx / 3 in this region. For / > 2z/ max the 
spectral slope is again 2/3. Reality will be some combination 
of these histories. 

We therefore expect the cosmological spectrum we cal- 
culate later (section 6) to be composed of a superposition 
of curves of these shapes, modified for chirp mass varia- 
tions, redshift effects and time delay between progenitor star 
formation and DD formation. The detailed calculations de- 
scribed in following sections follow in detail the evolution 
of all sources from ZAMS to merger, and do not rely upon 



approximate treatments of the kind given above. Simple es- 
timates of the background amplitude are discussed in section 
7. 



5 MODEL CONSTRUCTION 

5.1 The BSE code and population synthesis 

The rapid evolution code BSE (Hurley, Tout & Pols 2002) 
is used throughout this work whenever a binary system 
is evolved. This code is a fit to detailed models of stel- 
lar evolution, and produces an evolutionary time-sequence 
x(tj) of the properties x of any input ZAMS binary system. 
The code's time-resolution adapts to the shortest current 
timescale for change of the system components and orbit, 
due to e.g. nuclear evolution, angular momentum loss or 
mass transfer, which are all treated iteratively and have fi- 
nite duration. In this way, even the most fleeting of evolu- 
tionary phases is captured in detail, without requiring ex- 
cessive time resolution during long phases in which little 
changes. This is especially useful in the study of gravita- 
tional waves, since the majority of the GW emission from a 
given system occurs over an inspiral timescale much shorter 
than the nuclear timescales of the binary's parent ZAMS 
stars. Some of the most relevant features of the BSE code 
will be described in the following section; see Hurley et al. 
(2002) for full details. 

The output x(tj) from the code can be used to construct 
a stellar population at time T as follows. This method is 
similar to that used by Hurley et al. (2002) to characterise 
the Galactic binary population. 

We describe the ZAMS binary parameter space in terms 
of the primary (larger) mass Mi, secondary mass M2 (or 
mass ratio q — M2/M1 < 1), orbital semi-major axis a and 
orbital eccentricity e. We divide this space into grid boxes, 
and from each box k, we randomly choose a ZAMS system 
to represent the evolution of all sources in that box. 

The number Vk of sources born into box k per unit 
binary system realised is determined by probability distri- 
butions A(a), S(e) and $(Mi,M 2 ) = <E>i(Mi)/(q) in the 
ZAMS system properties described above (see section 5.3). 
Vk is obtained by integrating the product of these distribu- 
tion functions over the extent of box k. 

We wish to construct the population of sources present 
at time T. For each output timestep tj, the system with 
properties x(tj) can be viewed as a system born between 
times (T — tj+i) and (T — tj). If at this point the star for- 
mation rate was TZ — TZ(T — tj) (expressed as a number 
of binary systems born per unit time), then the number of 
systems with properties x(tj) we expect to see at time T is 
given by 

Nj,k{T) = {tj+i - tj)K(T - tj)V k , (12) 

so long as T > tj, so that stars were not born before time 
began. We perform this calculation for all boxes k and all 
timesteps j, so that the total population at time T is given 
by the combination of all Afj,k{T). 

This method of population synthesis ensures that 
sources from even unlikely regions of ZAMS parameter space 
are represented, weighted by their low formation probabil- 
ity. Coupled with the adaptive time-resolution of the BSE 
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output, and a sufficiently fine grid spacing, this technique 
allows the synthesis of a statistically reasonable population 
in a modest amount of computing time. Alternatively, statis- 
tical accuracy can be ensured with a Monte Carlo approach 
by simply generating a large enough number of stars under 
the initial distribution functions (see Belczynski et al. 2002). 

Our grid extends from 0.08 to 20 Mq in the mass of the 
primary Mi, and from 0.08 Mq to Mi in the secondary's 
mass. The initial separation is gridded from 2(Ri + R 2 ) to 
10 5 (7?i + R 2 ), where Ri and R 2 are the ZAMS radii of the 
primary and secondary respectively. We find that our back- 
ground fluxes are statistically accurate to around one per 
cent if we choose grid spacings of 0.05 in In M for each mass 
and 0.1 in In a for the separation. This corresponds to evolv- 
ing ~ 7 x 10 J binaries. For the Galactic tests described in 
section 5.3 we find that it is sufficient to use a grid spacing 
twice as large in each dimension. 

The BSE code has previously been tested against vari- 
ous Galactic populations of binary stars (Hurley 2000). A 
set of input parameters and distributions is recommended 
for use with the code, to best reproduce the observed Galac- 
tic binary population as a whole. However, in this work we 
are keen to quantify the effects of astrophysical uncertain- 
ties upon population synthesis calculations of the GW back- 
ground, and so in the following subsections we construct a 
set of models which differ in their choice of input parame- 
ters but produce specifically a Galactic DD population not 
in conflict with observations. The current observational un- 
certainties about DDs admit a range of models. This set of 
models is then considered representative of the population 
synthesis uncertainties affecting the GW background. 



5.2 The state of observations 

The observations of DD stars are currently undergoing a 
revolution. Full results of this revolution have not yet been 
published, so the detailed comparison of synthesised popu- 
lations with observations is still difficult. 

Marsh (2000) reported on the 15 then known DDs 
with measured periods, six of which had measured com- 
ponent mass ratios (Maxted, Marsh & Moran 2002). 
Searches for DDs have mainly focussed on low-mass WDs, 
A/wd<0.5 M (e.g. Marsh, Dhillon & Duck 1995), since 
these must have formed through giant stars losing their en- 
velopes in binary systems, before the helium burning that 
would inevitably occur in a single star. Maxted & Marsh 
(1999) determined that the fraction of DDs among these DA 
WDs is between 1.7 and 19 per cent, with 95 per cent con- 
fidence. Statistical comparisons with population synthesis 
models are thus difficult, given the sample size and level of 
bias, but there are some notable disagreements between ob- 
servations and theory that are not easily explained in terms 
of selection effects. The first of these is the lack of observed 
very low mass He WDs (M ~ 0.25 Mq). Theory predicts 
an abundance of such sources. Nelemans et al. (2001b) sug- 
gest that this can be explained by a more rapid cooling law 
for low-mass WDs than is commonly used. The second dis- 
crepancy is in the distribution of known DD mass ratios, 
which is seen to peak near unity (Maxted et al. 2002). Even 
considering selection effects (Nelemans et al. 2001b), this is 
difficult to explain in terms of either standard DD formation 



route, since as described in Section 3, the WD masses are 
expected to differ significantly. 

This prompted Nelemans et al. (2000) to suggest an 
alternative scenario in which a common envelope phase be- 
tween a giant and a main sequence star of similar mass does 
not result in a substantial spiral-in of the orbit, meaning 
that the second common envelope phase does not occur un- 
til the secondary's radius is larger (relative to that of the 
primary when it filled its own Roche lobe) than in the stan- 
dard CEE+CEE picture, so that the second WD formed is 
more massive, closer to the mass of the first-formed WD. 
They motivate this choice by parametrizing in terms of an 
angular momentum, rather than an energy balance (cf. sec- 
tion 3). 

The observational sample of DDs is currently being sub- 
stantially increased by the SPY project (Napiwotzki et al. 
2002), a spectroscopic study of ~ 1500 apparently single 
WDs (not restricted to low mass) to search for radial veloc- 
ity variations indicative of binarity. Napiwotzki et al. (2002) 
report that of the 558 WDs surveyed so far, 90 (16 per cent) 
show evidence for a close WD companion. Of these, mass 
ratio determinations are reported for three DDs (Karl et 
al. 2002), these three continuing the observed trend of mass 
ratios near unity. 

The results of the SPY project, once analysed fully, will 
help to constrain DD population synthesis calculations in a 
greatly improved way. However, given the preliminary and 
partial nature of the results so far, we can make only rather 
broad statements about their compatibility with any given 
synthesised Galactic population. This process is described 
in the next section. 

5.3 Candidate models 

Our fiducial population synthesis model (Model A) is similar 
to the preferred model suggested by Hurley (2000) (also his 
Model A): we use the initial mass function (IMF) of Kroupa, 
Tout & Gilmore (1993) (KTG) for <E>(Mi), we distribute M 2 
uniformly in the mass ratio q = M2/M1, f(q) — 1, and we 
start with a flat distribution in log a, choosing our limits 
as 2(Ri + R 2 ) < a < 10 5 (Ri + R 2 ), where Ri and R 2 are 
the ZAMS radii of the primary and secondary respectively. 
We have tidal effects switched 'on', we use a — 3.0 for the 
common envelope efficiency parameter, and we assign all 
stars solar metallicity, Z = 0.02. For the Galaxy, we adopt 
the constant star formation rate 1Z over the past 10 Gyr 
which gives a stellar disk mass of 6 x 10 10 M Q today. 

We differ from Hurley's Model A in three main ways: 
first, we assign an initial binary fraction of 50 per cent (cf. 
Hurley's 100 per cent) since this is observed locally to be 
the case (Duquennoy & Mayor 1991) and we evolve a set 
of single stars alongside the binaries, distributed according 
to the same IMF as the binary primaries. Second, we as- 
sign a ZAMS orbital eccentricity e to all systems, according 
to a thermal distribution H(e) = 2e,0 < e < 1.0. Hurley 
(2000) finds that an e = model gives a somewhat bet- 
ter fit to observations (though he finds that the numbers of 
close (P < 10 d) DD systems produced are not affected); 
we will also test a model of this type as part of our parame- 
ter variation (see below). Lastly, Hurley's Model A assumed 
the envelope binding energy parameter A = 0.5 for all stars, 
whereas here we allow this parameter to be calculated in 
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the code (values of A are from fits to detailed models of stel- 
lar evolution by O. Pols and are an addition to the code 
described in Hurley et al. 2002; J. Hurley, private communi- 
cation, 2003), and in addition we include 50 per cent of the 
envelope's ionisation energy in its binding energy. 

We test our synthesised Galactic populations against 
observations in a necessarily simple way. The aim is to re- 
ject models in clear conflict with the observed population 
of double degenerate stars, and to admit all others as rep- 
resentative of the uncertainties in DD population synthesis. 
Since the overall normalisation for the cosmological integral 
will be entirely separate from that used for the Galaxy, we 
choose primarily to compare relative populations as opposed 
to absolute numbers of Galactic sources. An ideal criterion 
is the fraction among field WDs of close DD binaries, which 
currently available SPY results place at 16 per cent. Since 
the sample size is substantially larger than that of Maxted & 
Marsh (1999), we adopt the SPY data, despite their incom- 
pleteness. We assume a negligible false-positive rate for SPY, 
and approximate the survey as magnitude-limited (V < 15) 
for the purposes of comparison. The somewhat approximate 
Galactic model and star formation history used here are suf- 
ficient, given the generosity of our selection criteria and the 
fact that we compare fractional quantities wherever possible. 

We distribute all stars according to a simple double ex- 
ponential Galactic disk model (scale height 200 pc, scale 
radius 2.5 kpc), then calculate the fraction of WDs with 
V < 15 expected to be members of DD binaries with 
P < 100 d. We then require that this calculated fraction 
be at least 10 per cent, if a given model is to be accepted. 
We assign a lower limit only, since our calculated binary 
fractions are likely to be overestimates, for several reasons. 
First, 100 d is a generous upper limit to the orbital peri- 
ods detectable with SPY; second, we do not address the 
issue of the substantial lack of observed low-mass (hence 
binary-member) WDs found in other population synthesis 
studies; and finally, the cooling curves used are the simple 
Mestel curves from Hurley et al. (2002); if we instead use 
the 'modified Mestel cooling' from Hurley & Shara (2003), 
which better fits the theoretical curves of Hansen (1999), 
then our calculated binary fraction decreases by a few per- 
cent. For our fiducial Model A, with Hurley et al. (2002) 
cooling, we find that 18 per cent of field WDs will show up 
as DDs in such a survey, in reasonable agreement with the 
SPY results. 

We also find a local total space density of WDs of 
9 x 10~ 3 pc~ 3 , and compare this with observational values, 
which range from ~ 4 — 20 x 10~ 3 pc~ 3 (Nelemans et al. 
2001b, and references therein). We do not attempt to com- 
pare to distributions in mass, mass ratio or period in detail: 
the observed distributions are subject to complex selection 
effects, and turn out often to be most constraining for WD 
cooling models (e.g. Nelemans et al. 2001b), whose devel- 
opment is beyond the scope of this paper. We note how- 
ever that in a volume-limited sense, the mean mass ratio 
(where q < 1 by definition) for detached WD-WD pairs is 
{q} = 0.62, not in good agreement with observations, but in 
common with other studies. 

We then go on to consider adjustments to our model, 
varying the initial distributions and mass transfer prescrip- 
tions. In all respects other than those mentioned below, 
these models are identical to Model A. 



In Models B, C and D, we use common envelope effi- 
ciency parameters a of 1.0, 2.0 and 4.0 respectively, while 
Model E uses the angular momentum formalism proposed 
by Nelemans et al. (2000) for the first phase of spiral in, 
with their recommended value of 7, and with a = 4.0. 

In models N, O, P and W, we also perturb the common 
envelope phase. In Model N, we include all of the envelope's 
ionisation energy (a positive quantity corresponding to the 
energy released when the ionised part of the envelope recom- 
bines) in its binding energy, meaning that envelopes will be 
less strongly bound and hence their removal will require less 
orbital shrinkage. This effect becomes important for stars 
on the AGB. Model O, on the other hand, does not include 
any of the ionisation energy. 

Determinations of A from stellar modelling are found 
to depend on the definition of the core-envelope boundary 
(Tauris & Dewi 2001) in giant stars. Because of this un- 
certainty, we also evolve models W and P in which we fix 
A = 0.5, with a = 3 and a — 4, respectively. 

In Model F, we choose the primary mass from the IMF 
of Scalo (1986), as in Schneider et al. (2001). Then in Model 
G we select both Mi and M2 independently from the KTG 
IMF, as suggested by Kroupa et al. (1993). We also evolve 
a Model K, in which initial orbital eccentricities are set to 
zero. 

Models L and M alter the production of DDs via the 
RLOF+CEE route described in section 3. It has been sug- 
gested (Han, Tout & Eggleton 2000) that Roche lobe over- 
flow may be stable until later in the Hertzsprung gap (HG) 
than happens using the BSE code, so a Model with enhanced 
HG overflow was added (Model L). Model M has semicon- 
servative overflow during this stage, to emphasise the uncer- 
tainties associated with HG mass transfer. 

The Galactic DD population was simulated using each 
model in turn; the results of this exercise are summarised 
in Table 1. Imposing the criterion given above, we eliminate 
Models B, G and W based on their under-production of DDs. 
If we increase the binary fraction to 100 per cent, this tends 
to under-produce single WDs, leading to an especially high 
DD fraction and a low overall WD space density. Note that 
the table also contains a Model H, which is in agreement 
with observations and is described in the next section. 

Thus the models A, C, D, E, F, H, K, L, M, N, O and 
P progress to the next round, as representative of reason- 
able astrophysical uncertainties in our population synthesis 
calculations. Three further models are added later (section 
6.4); these vary in their cosmic star formation and metallic- 
ity histories, and so cannot be tested against the Galactic 
DD population. 

5.4 Interacting DDs and modifications made to 
BSE code 

Some modifications were made to the BSE code regarding 
the treatment of accreting DD systems. In this we mainly 
follow the recommendations made in the detailed population 
synthesis work of Nelemans et al. (2001a). 

AM CVn stars are mass-transferring compact binaries 
in which the transfer is driven by gravitational radiation, 
and in which the accretor is a white dwarf and the donor is a 
Roche-lobe filling star, which could be another (less massive) 
white dwarf, or a helium star. For a review, see Nelemans 
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Table 1. Properties of Galactic DD models; details of models 
given in section 5.3. % DD is the percentage of field WDs in a 
magnitude-limited survey that will have a WD companion in an 
orbit with P < 100 d. pwD,0 is the local space density of WDs 
(single and double). {q) vol is the volume-limited average detached 
DD mass ratio q, where q < 1 by definition. 



Model 


% DD 


(~i~\~\TT\ r-\ 

(10- 3 pc- 3 ) 


(a) , 

W/vol 


Acceptable? 


A 


18 


9 


0.62 


Yes 


B 


7 


8 


0.68 


No 


C 


13 


9 


0.63 


Yes 


D 


20 


9 


0.63 


Yes 


E 


24 


9 


0.75 


Yes 


F 


22 


6 


0.64 


Yes 


G 


6 


6 


0.58 


No 


H 


18 


9 


0.62 


Yes 


K 


17 


9 


0.63 


Yes 


L 


18 


9 


0.63 


Yes 


M 


17 


9 


0.62 


Yes 


N 


17 


9 


0.63 


Yes 


O 


20 


9 


0.59 


Yes 


w 


9 


8 


0.62 


No 


p 


12 


8 


0.62 


Yes 



ct al. (2001a) and references therein. While not expected to 
be the dominant source of the Galactic gravitational wave 
background (Hils 1998; Hils & Bender 2000), some of these 
systems will be useful as 'verification' sources for LISA, with 
large, predictable gravitational wave amplitudes. 

We include in our definition of AM CVns all systems 
in which a helium star or WD is transferring mass on to a 
WD, including those systems in which the donor star is a 
CO or ONe WD. 

5.4.0.1 The WD family When the donor star is a white 
dwarf, the orbital separation at initial Roche lobe overflow 
is around 0.1 R,0, which is often sufficiently small that the 
accretion stream impacts directly on the accretor's surface, 
so an accretion disc is not expected to form. This has im- 
plications for the orbital evolution of the mass-transferring 
binary. When an accretion disc is present, tidal torques on 
the outer edge of the disc return to the orbit the angular 
momentum carried away from the donor by the accretion 
stream. In the absence of such a mechanism for restoring 
the orbital angular momentum, the criterion for stable mass 
transfer becomes much more stringent, and in most cases 
an AM CVn star will not form, precluding the existence 
of the WD family. Here we take the optimistic view (as in 
model II of Nelemans et al. 2001a) that, even if no disc 
is present, some tidal mechanism has an equivalent effect 
and that all WD-WD systems for which the mass ratio is 
< 0.628 (Hurley et al. 2002) will commence stable mass 
transfer upon Roche lobe overflow. We modify the BSE code 
accordingly. This optimism is perhaps warranted, since we 
do see WD family AM CVn systems, e.g. Israel et al. (2002) , 
which reports on the discovery of a helium-transferring com- 
pact binary with orbital period (321 s) too short to involve 
a (non-degenerate) helium star donor. 

5.4.0.2 The helium star family In this case, the donor 
star is a helium star, produced when a star with mass 



> 2 Mq loses its envelope on the RGB. Since these stars 
can live for a rather long time compared with the main se- 
quence lifetimes of their progenitors, there is a significant 
chance that through GW losses (or sometimes radial evo- 
lution) they will commence mass transfer before evolution 
to the WD stage. Here we shall employ the same condition 
on the dynamical stability of this mass transfer as Nele- 
mans et al. (2001a): q = M nH e/M wd < 1.2 (we use 'nHe' to 
denote (naked) helium star, to avoid confusion with helium- 
core WDs). Stellar modelling (Savonije, de Kool & van den 
Heuvel 1986) indicates that rapid mass transfer forces the 
helium star out of thermal equilibrium, increasing the ther- 
mal timescale beyond a Hubble time. The star cannot ever 
regain thermal equilibrium, and becomes semi-degenerate 
(as opposed to fully degenerate) as its mass falls. This re- 
sults in a negative exponent in the mass-radius relation, so 
that the orbital separation then increases as the helium star 
stably loses mass, i.e. an AM CVn system is formed. Note 
that at the onset of Roche lobe overflow, helium stars are 
always large enough that an accretion disk can form. 

The standard BSE code does not incorporate the pos- 
sibility of these semi-degenerate helium stars, so this was 
added. Here we adopt the same semi-degenerate mass-radius 
relation as in Nelemans et al. (2001a) (in solar units): 

i? nHc = 0.043 M-° c 062 , (13) 

and switch between this and the regular non-degenerate re- 
lation by selecting the larger of the two radii when the he- 
lium star is transferring mass on to a WD companion. In our 
code, this changeover occurs at M n u c ~ 0.29 Mq. We also 
modify the mass transfer rate prescription in the code, in 
order that the transfer responds more quickly to the initial 
overflow, so that the helium star does not hugely overhang 
its Roche lobe, and we halt further helium burning, so that 
the star cannot evolve to the WD stage during transfer, due 
to its long thermal timescale. We note that this modifica- 
tion is fairly crude, but ought to give a good indication of 
the relative importance of helium star AM CVn systems as 
sources of the GW background. 

A further issue in the formation of any helium- 
transferring system is that of edge-lit detonations (ELDs), 
which are believed to occur after a layer of helium has built 
up in the surface of an accreting CO WD. The BSE code det- 
onates CO WDs in this way after the accretion of 0.15 Mq 
of helium. We evolve separately a model (Model H) in which 
this is increased to 0.3 Mq, as in Model II of Nelemans et 
al. (2001a). 



6 COSMOLOGICAL EQUATIONS 

In this section we describe our calculation of the cosmolog- 
ical background. We adopt a standard lambda-cosmology, 
with n m = 0.3, n A = 0.7 and H = 70 km s" 1 Mpc" 1 . This 
means that the current age of the universe, To = 13.5 Gyr. 
We assume isotropy throughout; for an analysis of the small 
anisotropy due to the localisation of binary stars in galaxies 
which follow the large scale structure of the universe, see 
Kosenko & Postnov (2000). 
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6.1 Basic Equations 

The specific flux Ff r = dF(f T )/df T received at frequency f T 
from an object at redshift z with specific luminosity L/ e is 
given by (e.g. Peacock 1999) 



Ffr 



47Td L (z) 2 \df r 



(14) 



where f a — (l + z)f T , dh(z) = (l + z)divi(z) is the luminosity 
distance to redshift z and c?m is the proper motion distance 
(cf. section 5 of Hogg (2000), which is also l/(27r) times the 
proper ('comoving') circumference of the sphere about the 
source which passes through the earth today) . 

If the radiation comes from a large number of sources 
spread over redshift and isotropically distributed on the sky, 
we can write dLf e (z) = lf c (z)dV(z), where £f e (z) is the 
comoving specific luminosity density (say in erg s _1 Hz - 
Mpc -3 ), dV(z) = And^dx is the comoving volume element 
and x is the comoving distance. 

We can then write the specific flux received in gravita- 
tional waves today as 



J z = 



=0 4ndl(z) 

To 



tf.{T) fdf, 



=o (l + *CO) 



dV(z) 



c dT, 



(15) 



(16) 



using d\ = —(1 + z)cdT, where T is cosmic time. 

This is the basic equation on which the code is based. 
The equation is discretised in / r , T and £ as described in 
section 6.2. 



6.2 Computational Equations 

In the code, we bin the received gravitational waves in fre- 
quency. To calculate the flux received in a frequency bin 
with limits f T i and f r 2, we integrate Eq. 16 between these 
limits: 

r/r2 fT 

cdTdfr 



(T) fdf, 



+ z(T)) \df, 
To r(i+z)f t 2 £j (y") 

T= i(l+,)/ rl + 



df e cdT, 



(17) 



i.e. we integrate only over those emitted frequencies that will 
have been redshifted to arrive in this frequency bin today. 
The bin size was chosen to be 0.1 in log 10 (/ r ). 

Clearly, to calculate F, we need to know the comoving 
luminosity density £f e in gravitational radiation at frequency 
/ e as a function of cosmic time. 

We first obtain the source population at a given cosmic 
time Ti, by simply generalising Eq. 12, so that now 



N k ,j(Ti) = (t j+ i - tj)n c (Ti - tj)V k , 



(18) 



where 1Z C (T) is the cosmic star formation rate at time T, 
expressed as a number of binary stars born per unit time per 
unit volume, and Nk.j(T) is the number density of binaries 
with parameters k, j at cosmic time Ti, and where we require 
Ti > tj. 

The gravitational wave luminosity density at time T, is 
then given by 



(19) 



i.e. we simply sum over the emission at frequency f e from 
all sources k,j present at that time, weighted by their space 
densities. 

Since each binary source s emits radiation at only spe- 
cific frequencies f n — nv s (where v s is the orbital frequency 
of binary s) at a given time (Eq. 3), this sum can be ex- 
pressed as 

= E W w( T i) E 9(n, e Kj )5{f c -nv k ^).{2Q) 

k.j n 

We then have 

Ff^ fr2 = E E "E ^m^^h^ cAT, (21) 



i k,j n m i n 



(l + Zi) 



where we have also discretised the integral over cosmic time 
T, as a sum over i intervals AT, and where n is an integer, 
with the limits n m i n and n max defined by / r i < < / r2 - 

At a given redshift Zi{T), we just sum over those harmonics 
of those sources that will lead to emission at frequencies f c , 
with / r i(l + Zi) < f c < /r2(l + «i), and hence reception in 
the / r i — » / r 2 frequency bin today. 

The integration timestep AT must be sufficiently small 
that the emitting source population does not change sig- 
nificantly on timescales shorter than this, i.e. we assume a 
quasi-steady state population during this interval, so that 
our snapshot of the population at time T is representative 
of the whole timestep AT. A value of AT = To/50 was used 
throughout. We checked that timesteps smaller than this did 
not yield noticeably different results. Individual sources may 
evolve significantly within this timestep, but the character- 
istic emission of the population will be unchanged. It should 
also be noted that the evolutionary timesteps taken for the 
binary stars are independent of this integration timestep (see 
section 5.1), so that AT may be made much larger than the 
timescales of the evolutionary processes of interest, so long 
as the population is roughly steady-state over AT. 

Equation 21 is the sum performed by the code written 
for this paper, for a large number of received frequency bins 
over the range 10 -6 < f T < 10° Hz. For practical purposes, 
the sum over harmonics is truncated when g(n, e) drops be- 
low 10~ 3 , well beyond the peak in the emitted spectrum at 
n v — 1.63(1 — e)~ 3 / 2 . For typical e < 0.95, our numerical 
cutoff at g = 10~ 3 corresponds roughly to including only 
n < 5n p . The higher values n > 5n p contribute less than 1 
per cent of the total gravitational wave luminosity. 



6.3 Quantities Used 

Some quantities commonly used in gravitational wave 
astronomy are: T/ r (/ r ), 4tv times the specific intensity; 
f2 gw (/r), the fraction of closure density per logarithmic GW 
frequency interval; and the power spectral density Sh(fr)- 
The first of these, T/ r (/ r ), can be calculated from 



*>,(/«) = 



F /,l^f,: 



(22) 



(/rt-/rl)' 

The second, Q gw (f r ), is the fraction of closure energy 
density contained in gravitational waves received in the log- 
arithmic frequency interval around f T , i.e. 

frFf(fr) 



^gw(/r) — 



(23) 
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where p c is the critical mass density of the Universe; p c = 
(3Hi/8nG) ~ (1.88 x 10~ 29 ) h\ w g cm' 3 , where H = 
100/noo km s _1 Mpc -1 . In terms of computational quan- 
tities, 

~ 0.0175 F ftl ^ ft2 (24) 

(where -F/ rl ^/ r2 is in erg s _1 cm~ 2 ) since A(ln/ r ) = 
ln(10° - 1 ) ~ 0.23 and h wo = 0.7. 

The power spectral density Sh(fr) is given by 

Sm = ^Ffr(fr). (25) 

f 1/2 

Usually this is plotted as S\ /2 ~ (5.6 x HT^-^-Hz- 1 / 2 , 
where Ff r is in erg s _1 cm -2 Hz -1 , and / r is in Hz. 

6.4 Cosmic star formation history 

As pointed out by Schneider et al. (2001), most determina- 
tions of cosmic star formation history are based on the UV 
emission from massive stars (e.g. Madau et al. 1996; Steidel 
et al. 1999), and use an assumed single-star IMF (commonly 
that of Salpeter 1955) to convert observed UV flux into a 
star formation rate as a function of redshift. This type of 
rate is inconvenient here for two reasons: first, a non-trivial 
factor is required for conversion to a binary star formation 
rate (for an assumed binary fraction), because of the need to 
correct for the observed flux from companion stars; and sec- 
ond, the total star formation rate is pivoted on the high-mass 
end of the stellar distribution, while here we are interested 
in studying the remnants of low- to intermediate-mass stars. 
This results in a crucial dependence on the choice of stellar 
IMF. 

Schneider et al. (2001) overcome the first problem by 
assuming the measured shape of the cosmic SFH as a func- 
tion of time, but normalising its amplitude to the local rate 
of core-collapse supernovae. This Type Ibc/II SN rate is a 
more easily calculated quantity for a given (binary or sin- 
gle) IMF than is the UV luminosity density. Since Schnei- 
der et al. (2001) are also concerned with neutron stars in 
their study, this is a reasonable choice. However, the second 
problem remains when one is concerned with WDs; and in 
addition, not only does the normalisation pivot on the high 
mass stars, but it also depends crucially on the ratio of lo- 
cal to peak cosmic SFR. We also note that the minimum 
mass of star producing a core-collapse supernova explosion 
is uncertain (e.g. Jeffries 1997). 

For our normalisation, we use instead the observed local 
stellar mass density Q t , as derived from the local near-IR 
luminosity function by Cole et al. (2001). This quantity is 
most sensitive to stellar masses near the MS turnoff in old 
populations, M ~ 0.8 — 1.0 Mq, and thus is more closely re- 
lated than the SNIbc/II rate to the DD progenitor popula- 
tion. We convert between their assumed single star IMF and 
our binary star IMFs by keeping constant the mass in stars 
in this range. We then use the recycled fraction R = 0.42, 
as for the Kennicutt (1983) IMF used in Cole et al. (2001), 
to convert stellar density today to total mass of stars ever 
formed, Q*,tot (the time- integral of the cosmic star formation 
rate). Doing this, we obtain Q*,tot = 5.0 x 10 8 M Mpc~ 3 
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z 

Figure 3. Two possible cosmic star formation histories, plotted 
as a function of redshift z and parametrized according to the 
smooth curve fits given in Cole et al. (2001). Dashed line: no 
extinction correction made, used in Model J. Solid line: extinction 
corrected with E(B — V) = 0.15, used in all other Models. The 
time integral of each rate is fixed using the appropriate f2*,tot 
derived from Cole et al. (2001). Curves shown are for the KTG 
IMF. 



for the KTG IMF, while for the Scalo IMF this figure is 
fl*,tot = 4.0 x 10 s M Mpc~ 3 . Due to this rather crude con- 
version, the uncertainty in these figures will be greater than 
the 15 per cent quoted by Cole et al. (2001) for Q„; we esti- 
mate the resulting uncertainty to be ~ 30 per cent. 

Cole et al. (2001) note that their calculated stellar den- 
sities are most consistent with UV-derived star formation 
rates if the extinction corrections used in these methods are 
moderate. However, we would like to assess the effects of un- 
certainty in the shape of the cosmic star formation history. 
We therefore select both a history with large extinction cor- 
rections and one with none, keeping the integral over time 
fixed to f2* jto t for each. The corresponding curves are plot- 
ted in Fig. 3. We use the extinction-corrected rate, favoured 
by Steidel et al. (1999), in Model A and all other models 
except for Model J, which uses the uncorrected rate (but is 
identical to Model A in all other respects). We also intro- 
duce Models Q and R, whose metallicity histories differ from 
that of Model A: in Model Q, stars born during the first Gyr 
have metallicity 1 /20 solar, while stars born later have solar 
composition; in Model R, all stars have metallicity Z = 0.01, 
i.e. half-solar. 



7 BASIC RESULTS 

The GW background spectrum received in the frequency 
range 10~ 6 < / r < 10' 1 Hz, generated using our fiducial 
Model A, is plotted in Fig. 4. The total amplitude is bro- 
ken down into separate contributions from four main evo- 
lutionary stages: main sequence-main sequence (MS-MS), 
WD-MS, WD-WD and WD-helium star (WD-nHe) bina- 

1/2 

ries, and plotted in terms of each of Ff T , Q gw and S h de- 
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Figure 4. The GW background for our fiducial Model A, in terms 
of the three quantities described in section 6.3. Solid line: WD- 
WD pairs; dotted line: nHe-WD pairs; dashed line: MS-MS bi- 
naries, and dot-dash line: WD— MS binaries. The total signal (the 
sum of the four parts) is given by the thick solid line. Only n = 2 
harmonics of the orbital frequency are plotted (see section 7.1.1). 

scribed in section 6.3. The unitless O gw will be our preferred 
quantity for the remainder of the paper 1 . 

1 Note that since Cole et al. (2001) quote f2*/i in their paper, 
and we use this quantity to normalise our star formation rate, 
our calculated f2 gw also scales as h —1 . We use h = 0.7.) 



10 " 




f'Hz 

Figure 5. Comparison of spectral shapes for all Models. Curves 
are the same as in Fig. 4, but all quantities are plotted as solid 
lines. Only n = 2 harmonics of the orbital frequency are plotted 
(see section 7.1.1). 



The four component spectra are plotted in Fig. 5 for all 
of the models evolved, to illustrate that the spectral shapes 
are largely unaffected by any of the changes made. A sum- 
mary of important quantities for each model (to be discussed 
later) is given in Table 2. For reference, we also list a Model 
A', identical in parameters to Model A, to demonstrate the 
typical level of statistical variation in the results. This is 
clearly at the 1 per cent level in flux, so that variations larger 
than this between models can be ascribed to parameter, and 
not statistical, variations. 

Throughout we will focus on the properties of the spec- 
trum around 1 mHz, in the centre of the LISA band and 
of the spiral-in regime. We will also compare with the spec- 
tral properties at 10 mHz, at which frequency lower-mass 
WD- WD pairs can no longer be present and at which point 
this extragalactic WD-WD background will be the domi- 
nant LISA background source (see Fig. 16). 

It is clear that the signal in the LISA frequency band 
(0.1 < fr < 10 mHz) is dominated by the WD-WD compo- 
nent, as expected. Neither the MS-MS nor the MS- WD bi- 
naries can radiate at frequencies above the bottom of this 
bandpass, since even the lowest mass MS stars come into 
contact at frequencies below 1 mHz. WD-nHe pairs can con- 
tribute to a somewhat higher frequency due to the smaller 
radii of helium stars, but still come into Roche lobe contact 
at / e ~ 1 mHz. 

The WD-WD component clearly displays the spectral 
shape predicted in Section 4 (f2 gw oc / r F/ r , plotted in Figs. 1 
and 2), with a clear separation between static and spiral- in 
regimes at around 10 -4 Hz. The slope in the static regime 
suggests that sources are injected with a spectrum closer 
to iVb(^) oc v~ 2 / 3 than to but agreement to this level 
is encouraging. The spiral-in slope is slightly steeper than 
predicted, but this is due to the spectrum seen being the sum 
of spectra from populations with different chirp masses, as 
well as different merger and maximum injection frequencies 
(see Fig. 8), whose individual slopes in the spiral- in regime 
are closer to the predicted 2/3. Agreement with our simple 
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Table 2. Summary table of results from all models described in 
the text. f2g w (l mHz) is in units of 10~ 12 . R+C refers to the 
RLOF+CEE route to the DD stage, while C+C refers to the 
CEE+CEE route (see section 3; note that for Model E, we hold 
the RLOF+CEE contribution fixed from Model A). The flux- 
weighted mean chirp mass (M) contributing at f r ~ 1 mHz is in 
units of Mq , the next column lists the percentage contribution to 
O gw at 1 mHz from interacting binaries and the last column gives 
the inspiral remnant density No today, in units of 10 6 Mpc — 3 . 
Models B and G were rejected for reasons noted in Table 1. 



Model 


n 


gw (l mHz) 


(M) 


% 


N 




Total 


R+C 


C+C 




AM CVn 




A 


3.57 


1.35 


2.22 


0.45 


13 


1.17 


A' 


3.61 


1.36 


2.26 


0.45 


14 


1.18 


C 


3.06 


0.60 


2.47 


0.44 


16 


0.90 


D 


3.66 


1.64 


2.02 


0.47 


13 


1.20 


E 


4.21 


1.35 


2.86 


0.47 


10 


1.57 


F 


1.94 


0.72 


1.22 


0.41 


13 


0.75 


H 


4.10 


1.53 


2.58 


0.43 


25 


1.17 


J 


3.62 


1.38 


2.24 


0.45 


13 


1.17 


K 


4.29 


2.09 


2.20 


0.48 


13 


1.29 


L 


3.80 


1.53 


2.27 


0.45 


12 


1.25 


M 


2.80 


0.66 


2.14 


0.44 


15 


0.92 


N 


3.43 


1.36 


2.07 


0.46 


13 


1.13 


O 


3.89 


1.31 


2.57 


0.46 


13 


1.27 


P 


5.46 


1.00 


4.46 


0.55 


16 


1.20 


Q 


3.73 


1.43 


2.30 


0.44 


13 


1.32 


R 


3.83 


1.48 


2.35 


0.44 


12 


1.28 



predictions is therefore good and we feel that we understand 
well the origins of the spectrum. 

7.1 Contributors 

The breakdown of contributions to the background received 
at 1 mHz for our fiducial Model A is given in Table 3. In 
this section we identify the dominant source types, and those 
types whose contribution is negligible, then attempt to char- 
acterise the emitting population in terms of a mean chirp 
mass and inspiral remnant density. 

7.1.1 Eccentric harmonics 

As described in section 2, systems with eccentric orbits emit 
gravitational waves at all harmonics nv of the orbital fre- 
quency, not just the n = 2 harmonic as for circular orbits. 

The only close binaries we expect to be eccentric are 
unevolved MS-MS binaries in which tidal forces have not 
yet circularised the orbit. Almost every close evolved (e.g. 
WD-MS, WD-WD) system will have at some point expe- 
rienced a Roche lobe-filling phase, which will likely have 
circularised the system, through tidal circularisation and/or 
common envelope evolution. Figure 6 shows the contribution 
from harmonics with n > 3 to the MS-MS GW spectrum for 
Model A (which has a thermal initial eccentricity distribu- 
tion). Clearly the n > 3 harmonics contribute < 10 per cent 
of the MS-MS spectrum at frequencies / r < 0.5 mHz, and 
although they dominate the MS-MS spectrum above this 
frequency, these signals are buried deep below the other con- 
tributors at / r > 0.5 mHz (see Fig. 4). Hereafter we safely 
neglect the n 2 contributions to L gw , in the interests of 



Table 3. Percentage contribution to Q gw at 1 mHz from different 
DD pairs, for both contribution to total integrated background 
and contribution to background coming from the local universe, 
z = 0. All for fiducial Model A. MS-MS and WD-MS binaries 
contribute negligibly at this frequency. ('nHc' denotes a naked 
nondegenerate helium star.) All contributions, including the AM 
CVn values, are given as fractions of the total flux at 1 mHz. 



Pairing 


% over all time 


% locally 


He-He 


12.4 


29.5 


He-CO 


23.0 


25.3 


He-ONe 


0.6 


0.6 


CO-CO 


42.2 


33.2 


CO-ONc 


8.1 


4.4 


ONe-ONe 


1.0 


0.2 


(of which AM CVn) 


3.6 


4.7 


nHe-WD 


12.7 


6.9 


(of which AM CVn) 


9.7 


2.0 


Total 


100 


100 


(of which AM CVn) 


13.3 


6.7 




10 5 10 4 10 3 



f / Hz 

Figure 6. The GW background from harmonics with n > 3 from 
MS— MS pairs (thin solid line), plotted along with the total MS- 
MS pair contribution (dashed line), and the total background 
from all sources (thick solid line), demonstrating that the har- 
monic contribution is negligible. All for Model A. 

computing time, though we do not neglect eccentric orbits 
in computing stellar evolution sequences. 

7.1.2 Interacting binaries 

Interacting binaries (those in which either a WD or nonde- 
generate naked helium (nHe) star is transferring mass on to 
a WD) contribute 13 per cent of the GW background at 1 
mHz in Model A. Since at this frequency the majority of nHe 
star companions fill their Roche lobes, most of the nHe-WD 
background comes from interacting systems. At 10 mHz, 26 
per cent of the GW signal comes from interacting binaries, 
all of these necessarily WD-donor systems. The GW spec- 
trum due to interacting binaries is compared with the total 
signal in Fig. 7. 

The percentage contribution from interacting systems is 
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fairly constant across models, except for Model H, in which 
an accreting CO WD is permitted to accumulate 0.3 M Q 
of helium before detonation, as opposed to the 0.15 Mq in 
our fiducial model. This increase in survival rate boosts the 
interacting binary signal at 1 mHz by a factor of two. For 
the other models, the interacting WD-WD signal is boosted 
when the WD-WD pairs formed typically have larger mass 
ratios, so that more systems can commence stable transfer 
upon Roche contact, e.g. Model C. 

The spectral shapes from interacting systems are gov- 
erned by the mass-radius relation of the Roche lobe-filling 
star, and so do not share the spectral slopes displayed by the 
detached binaries. The overall contribution from interacting 
pairs is sufficiently small, however, that the total spectral 
shape is little affected by their presence. This is in line with 
results for the Galaxy found by Hils & Bender (2000) and 
Nelemans et al. (2001c). 

We can predict the spectral shape due to interacting 
WD-WD binaries using some simple scaling relations (in the 
notation of section 4): for a Roche lobe- filling WD of mass 
M d , we have M~ 1/3 oc R d = R L oc aM d 1/3 oc M\' 3 f' 2 ' 3 , 
using Kepler's law (for conservative mass transfer). If we 
then assume that the mass of the donor WD is much less 
than that of the accretor, then the system chirp mass M oc 
M^ 5 , so that / oc _M 3//5 and the system gravitational wave 
luminosity L gw oc f 10 / 3 M 10/3 oc f 16 ' 3 . 

For sources sweeping (backwards) through frequency 
space, we have N(f) oc 1// oc M~ 6/3 f- 11/3 oc /" 14/3 . 

Putting these together, we then have, for the emit- 
ted flux in the logarithmic frequency interval around /, 
figw(/) oc fF(f) = /L gw iV(/) « f' 3 . From Fig. 7, we 
measure the spectral slope between 0.4 and 6 mHz to be 
~1.7, in good agreement with this calculation. Interacting 
WD-WD sources are not present below this frequency range 
because evolution to these frequencies requires more than a 
Hubble time. Above ~ 6 mHz, the spectral shape depends 
on the fraction of sources of high enough mass to radiate at 
a given frequency; this number drops rapidly with increas- 
ing frequency. Note that, within the 0.4-6 mHz range, since 
the spectrum O gw oc / 5/ ^ 3 for interacting WD-WD binaries 
rises relative to f2 gw oc f 2 ^ 3 for inspiralling detached bina- 
ries, interacting binaries are more important contributors at 
high frequencies than at low. 



7.1.3 WD types, chirp mass and merger rates 

The dominant component of the background at frequencies 
0.1-10 mHz comes from the inspiral of WD-WD systems. 
From Table 3, we see that approximately half of this back- 
ground comes from CO-CO pairs, descended primarily from 
higher mass progenitors than the majority of He-He sys- 
tems. The dominance of these systems is a result of both 
the shorter time delay between star formation and DD birth 
for more massive MS stars, and the larger chirp masses for 
CO-CO systems, since the flux in the inspiral part of the 
spectrum scales as fi gw oc f 2 ^ 3 M 5 ^ 3 (see section 4). These 
two factors outweigh the fact that, from the IMF, many 
more potential progenitors of He WDs are born than those 
which always produce CO or ONe WDs after envelope loss. 

Figure 9 shows however that, as more low-mass MS 
stars evolve to the DD stage, the relative contribution to 
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Figure 7. The spectrum due to interacting binaries, for our fidu- 
cial Model A. The solid line shows the total WD-WD binary 
contribution, while the dotted line gives the spectrum from in- 
teracting binary WD-WD systems. The dashed line is the total 
nHe-WD spectrum, of which the dash-dot line gives the interact- 
ing nHe-WD contribution. 
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Figure 8. The background received from different WD-WD pair- 
ings, for Model A. From top to bottom at 1 mHz: thick solid line: 
CO— CO, thin dashed line: He-CO, thin solid line: He-He, thick 
dashed line: CO-ONc, thick dotted line: ONc-ONc, thin dotted 
line: He-ONe. 



the GW luminosity density from pairs involving He WDs is 
rising, and will eventually dominate. The percentage con- 
tribution to the local (z = 0) WD-WD GW emission at 
1 mHz from pairs including at least one He WD is 55 per 
cent, whereas their contribution to the integrated cosmolog- 
ical background received today is only 36 per cent. 

A useful way to look at this is through the chirp mass 
distribution. Shown in Fig. 10 is the contribution to f2 gw 
at 1 mHz as a function of system chirp mass (defined in 
section 4) for Model A, giving a flux- weighted mean chirp 
mass of 0.45 Mq. As increasingly lower mass systems evolve 
off the main sequence and become close DD pairs, this mean 
chirp mass is decreasing with time, as shown in Fig. 12. The 
chirp mass distribution depends on GW frequency (Fig. 11), 
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Figure 9. The contribution to f2 gw (l mHz) received today as 
emitted from each shell of cosmic time, AT = To/50, from each 
source type, for Model A. Linestylcs arc as in Fig. 8, with the 
addition of the thin dash-dot line for nHc-WD pairs. 
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Figure 10. Relative contribution to Q gw at 1 mHz as a function 
of chirp mass, for Model A, giving a mean flux-weighted chirp 
mass (M) = 0.45 M . 

most notably shifting towards higher masses at frequencies 
above which lower-mass WD-WD pairs will have merged. 
The mean chirp mass is somewhat higher below the critical 
spiral- in frequency, since for f c < 2f cr it, we have O gw oc 
M 10/3 , and above 2f cr i t , Ogw oc M°/ z (see section 4). 

Phinney (2002) derived a simple expression for the GW 
background in terms of the chirp mass M, assumed con- 
stant across all sources, and the current space density iVo 
of remnant spiralled-in sources (with a weak dependence on 
cosmology and star formation history). We can assess the 
usefulness of this formula as a predictor of the background 
flux by using the results of our population synthesis calcu- 
lations to see whether the computed fluxes can indeed be 
described by these two parameters only. 

To calculate the remnant density No, we first calculate 
the source spiral-in rate as a function of cosmic time. The 
rate of occurrence of Roche lobe contact between WD-WD 
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Figure 11. Flux-weighted chirp mass contributing to the GW 
background received, as a function of frequency, for Model A. 
Solid line: detached WD-WD pairs only, dashed line: all source 
types. The dip seen in this curve around 0.1 mHz is due to low- 
mass main sequence stars. 
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Figure 12. The flux-weighted mean chirp mass contributing to 
emission received today at 1-3 mHz (solid line) and 3-10 mHz 
(dashed line), from each shell of cosmic time. All for Model A. 

pairs (we shall call this the spiral-in rate) is different from 
the rate of WD-WD mergers, since for some subset of sys- 
tems (those with mass ratios q < 0.628) stable mass transfer 
will commence upon overflow, and an AM CVn system will 
form. We keep track of both of these rates here. 

For / greater than both 2i/ max and 2^ cr it, i.e. in the 
part of the spiral-in regime above which sources are born 
(see section 4), then for a quasi-constant spiral- in rate N 
over the timestep To/50, the continuity equation (Eq. 6) 
simplifies to 

N = ^vN h (26) 

i 

summed over all sources i at any given frequency satisfying 
the above requirement. For each source, v is given by Eq. 
7. We perform this sum at each step in cosmic time, using 
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systems with orbital frequencies in the range 0.8 < v < 
1.6 mHz, which is above the maximum injection frequency 
for the majority of sources, and below those frequencies at 
which the lowest mass WDs are coming into contact. We 
note that the inspiral time from v ~ 0.5 mHz is less than 
To/50 = 0.27 Gyr for all 7W>0.05 M Q , so that at each 
timestep we are accurately representing the spiral-in rate 
at that time. The only exceptions are very low chirp mass 
systems, which we neglect here anyway, since these will be 
interacting binaries, which are spiralling out. We also neglect 
all nHe-WD pairs, since the evolution of these systems is not 
governed exclusively by gravitational radiation, but also via 
radial evolution of the nHe star, and also because Roche 
lobe contact occurs for these systems within our frequency 
range. 

The spiral-in and merger rates obtained from Model A 
are plotted in Fig. 13. The present-day remnant density iVo 
needed for the formula of Phinney (2002) is the time integral 
of the spiral-in rate, since this gives the total number of 
sources that have contributed to the background. From our 
calculated rate, we obtain N = 1.17 x 10 6 Mpc -3 . 

Phinney (2002) deals only with the GW emission from 
non-interacting WD-WD systems, and so we should com- 
pare its predictions with only the non-interacting component 
of our computed signals, in addition to using a characteris- 
tic chirp mass M! for just those systems. For Model A, our 
flux-weighted mean chirp mass for detached WD-WD pairs 
is (M'} = 0.47 M at 1 mHz. Eq. 16 of Phinney (2002), 

converting to hioo = 0.7, and omitting the ^(1 + z) -1 ^ 3 ^) 
scaling factor in the interests of simplicity, becomes 

Using {M') and No for Model A in the above, we find 
f2 gw (l mHz) = 3.7 x 10~ 12 . We compare this with the com- 
puted value for detached WD-WD pairs, Sl gw (l mHz) = 
3.0 x 10~ 12 , and note that these agree to within 25 per cent. 
If we perform this same calculation for the other Models, we 
find that Eq. 27 overestimates the computed background by 
a similar fraction. 

The variation between models is thus well fitted by the 
formula. The relative fluxes are reproduced by Eq. 27 to 
within 5 per cent for all Models except D and E, whose fluxes 
relative to Model A are overestimated by 7 and 16 per cent 
respectively. The dominant scaling is due to variations in 
iVo, since in most cases {M') varies little between Models. 
For the cases in which {M') does significantly change (D, E, 
F, K and P), the omission of the chirp mass scaling in Eq. 27 
can improve (D, E) or worsen (F, K, P) the agreement with 
the results of our detailed calculations. This is perhaps as ex- 
pected, since our flux- weighted chirp mass is in fact not the 
same average as that required in the generalisation of Phin- 
ney (2002) to accommodate a range of chirp masses. Such a 
value would also incorporate the redshift-scaling omitted in 
the above. We note, however, that neither iVo nor either def- 
inition of M! is a directly observable quantity, requiring as 
they do integrations over cosmic time, and so are not easily 
determined from observations. 

The computed spectral shape is not precisely fi gw oc 
jV A (see Fig. 4), so we do not expect an exact reproduction 
of the spectrum using this formula. However, we conclude 
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Figure 13. The rate of WD-WD spiral-in as a function of cos- 
mic time. The thick solid line gives the total spiral-in rate, while 
the thin solid line shows the merger rate, that is the inspiralling 
sources that will merge, and not commence stable mass transfer 
(i.e. become AM CVn binaries) upon Roche lobe overflow. The 
thin dashed line gives the rate of merger of WD-WD pairs with 
combined mass > 1.4 Mq. For reference, the cosmic star forma- 
tion rate, multiplied by 1/(1000Mq), is plotted as the thin dotted 
line. All for Model A. 

that with a knowledge of No and M' , we can quickly predict 
the detached WD-WD background amplitude and to some 
extent its variation if these values change. We note however 
that a full population synthesis calculation enables the in- 
clusion of interacting systems, as well as the extraction of 
detailed spectral shapes and source property distributions, 
which are not available in a quick 'manual' calculation. 

7.2 Progenitors 

Here we outline the relative contributions from the two main 
pathways to the DD stage, and we assess the impact upon 
each of these routes of varying the population synthesis 
model. 

Figure 14 shows the contribution to fi gw at 1 mHz as a 
function of the initial mass of the primary, for Model A. The 
descendants of primaries with ZAMS masses in the range 2- 
4 Mq contribute 50 per cent of the signal, the flux-weighted 
mean progenitor primary mass being 3.7 Mq. Most of the 
sources in this range are the progenitors of CO WDs, since 
for M > 2 Mq , a CO WD will be produced via a helium star 
upon envelope loss on the RGB, and a CO WD will be pro- 
duced directly if the envelope is lost on the AGB. At 10 mHz, 
the mean progenitor mass rises to 4.7 Mq, since the (neces- 
sarily more massive) WD-WD pairs contributing there are 
descended from only the more massive ZAMS systems. The 
equivalent secondary mass distribution is not plotted here, 
but is always peaked towards initial mass ratios of unity. 

Of perhaps more interest is the distribution in initial 
orbital semimajor axis (Fig. 15, for Model A), which has a 
clear bimodal form, the peak at a ~ 5a m i n corresponding 
to DDs which formed via RLOF+CEE, and the peak at 

Qj ~ 50flrnin 

corresponding to the CEE+CEE route. We can 
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Figure 14. Contribution to f2 gw (l mHz) received, as a function 
of ZAMS mass of the primary, for our fiducial Model A. 




1 10 100 1000 10000 

a/2(R ]+ R 2 ) 



Figure 15. Contribution to f2 gw (l mHz) received, as a function 
of initial progenitor scmimajor axis, expressed as a ratio of the 
initial scmimajor axis to the minimum separation permitted in 
the code. For Model A. 

therefore approximately determine the relative contributions 
from these two routes by dividing this distribution between 
the two peaks (at a ~ 10a m i n for most models); the result 
of this division for each model is shown in Table 2. We note 
that the location of the CEE+CEE peak at a ~ 50a m i n and 
the typical masses of the dominant progenitor stars mean 
that for this route, the dominant pathway involves primary 
overflow on the AGB, followed by secondary overflow on its 
RGB. 

For Model A, fi gw (l mHz) = 1.4 x HT 12 (~ 38 per 
cent of the total) comes from sources that evolved via the 
RLOF+CEE pathway. Since the WD-WD pairs from this 
route are generally more massive than CEE+CEE pairs, the 
percentage contribution at 10 mHz from this route rises to 
44 per cent. 

In general, we shall find that it is the RLOF+CEE con- 
tribution that is affected more by varying the population 



synthesis model. Although it can be affected significantly 
by varying the form of the common envelope prescription 
(Models E and P), the CEE+CEE signal is quite robust to 
changes in the common envelope efficiency, since if systems 
originating at one separation happen to coalesce in a com- 
mon envelope phase, using a given model, there exists a shell 
of sources at greater a to take their place as the closest WD- 
WD systems at birth, out to a maximum of a ~ 10 3 o m i n at 
which Roche lobe overflow no longer occurs on the RGB or 
AGB. Webbink & Han (1998) describe this effect in terms of 
shifting the 'window' in initial parameter space from which 
the closest DD systems are descended. The weak dependence 
of results upon the common envelope efficiency parameter is 
also seen in population synthesis calculations for other types 
of binary, e.g. Kalogera & Webbink (1998) for LMXBs. 

Returning to the DD case, the RLOF+CEE pathway 
has no similar resource, occurring only in the rather narrow 
range of initial separation in which RLOF commences in the 
Hertzsprung gap. If we destroy more of these sources in the 
ensuing CEE phase, we lose more of the contributions from 
the RLOF+CEE route. 

Decreasing a (Model C) has this kind of deleterious ef- 
fect upon the RLOF+CEE pathway, but slightly increases 
the signal from CEE+CEE sources, since the systems which 
survive to the close DD stage were on average more widely 
spaced than for a = 3.0, so that the giant stars were physi- 
cally larger, i.e. more evolved, on average upon Roche lobe 
overflow, so gave rise to more massive WDs (also with more 
widely differing masses). This corresponds to moving the 
second peak in Fig. 15 to higher a. The lower mean chirp 
mass is largely attributable to the increased number of low- 
chirp mass interacting binaries present at this frequency, 
since the typical WD-WD mass ratio is larger, as described 
in section 7.1.2. Increasing the efficiency parameter (Model 
D) has the opposite effect upon each route. If on the other 
hand we use the common envelope formalism of Nelemans et 
al. (2000) (Model E), it becomes less simple to disentangle 
the two routes, since now they overlap somewhat in initial 
a— space, but since we know that this modification ought 
not to affect the RLOF+CEE contribution, we hold this 
fixed from Model A. The new CEE+CEE value turns out 
to be significantly enhanced, since a wider range of initial 
separations has been opened up to double common envelope 
survival. The nearer (by design) equality of WD pair masses 
leads to a decrease in the number of WD-WD AM CVn 
systems produced, and hence a smaller contribution from 
interacting systems than for Model A. 

The envelope ionisation energy becomes a significant 
part of the energy balance in AGB stars, and so its in- 
clusion is important in common envelope phases that com- 
mence at large orbital separations. Increasing the fraction of 
this energy included in the envelope binding energy (Model 
N) therefore decreases the number of wide binaries able to 
shrink enough to form close WD-WD pairs. Omitting it en- 
tirely (Model O), thus increasing the envelope binding en- 
ergy, has the opposite effect. 

Model P shows the greatest departure from Model A 
in terms of GW flux and mean chirp mass. The progenitor 
mass distribution for Model P is peaked towards higher mass 
(6-8 Mq) stars than for other Models. These differences 
can be traced to the outcome of common envelope phases 
on the Hertzsprung Gap (HG). The BSE fitting formula re- 
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turns values of A substantially smaller than 0.5 for most HG 
stars, corresponding to a high degree of central concentra- 
tion. Therefore using A = 0.5 in Model P results in much 
less shrinkage in these situations. 

High mass stars expand in radius by a large factor in 
their HGs, so that the final Roche contact (for both path- 
ways) is often a common envelope phase involving a HG 
star. The survival rate from this CE phase is boosted by the 
fixed lambda as described above. The resulting GW flux is 
therefore also greatly boosted for these higher mass stars, 
whose descendent WDs are sufficiently massive that rela- 
tively few are required to dominate the background GW 
flux. Given however that small values of lambda are robust 
for HG stars (they are also seen in the calculations of Dewi 
& Tauris 2000), we choose not to consider this prescription 
as a reasonable uncertainty on the background. 

The lower chirp mass (M.) seen for Model H is due to the 
inclusion of an increased number of interacting sources at 1 
mHz, compared with Model A; the value {M') appropriate 
for just detached WD-WD pairs for this model (used in the 
previous section) is the same as for Model A. 

Starting all systems with circular orbits (Model K) 
boosts the RLOF+CEE pathway, because fewer systems 
given initially tight orbits are lost due to immediate col- 
lision at periastron. Since systems descended from the 
RLOF+CEE route are generally higher-mass, the mean 
chirp mass for Model K is higher than for Model A. 
CEE+CEE route systems are little affected; the high-a peak 
in Fig. 15 is simply narrowed in a-space, since orbital sepa- 
rations are no longer altered by tidal circularisation before 
Roche contact. 

Aside from orbital circularisation, the main role of tides 
in the evolution to the DD stage is in orbital shrinkage before 
Roche contact, due to spin- up of the giant star. Neglecting 
tidal effects is thus similar to increasing the common enve- 
lope efficiency parameter, i.e. the progenitors of close DDs 
from the CEE+CEE route have smaller initial orbital sep- 
arations, and the DDs produced have smaller chirp masses 
on average. If on the other hand tidal effects were much 
stronger than in the BSE code, then we would expect little 
impact upon this route, since giant-star corotation is already 
typically achieved before Roche lobe overflow with the tides 
in BSE. 

The CEE+CEE route is as expected largely unaf- 
fected when we perturb dynamically stable mass trans- 
fer on the Hertzsprung gap. Much as one might expect, 
the RLOF+CEE route is enhanced when one enhances the 
transfer on the Hertzsprung gap (Model L), so that more 
mass is transferred to the companion, and more systems 
avoid a common envelope phase during the first phase of 
mass transfer (which tends to lead to merger). The orbit 
is also widened to a greater extent during transfer, mean- 
ing that more systems will survive the common envelope 
phase when the secondary evolves. Making the transfer semi- 
conservative (Model M) has an opposing effect; the orbit is 
widened less during stable overflow, meaning that more sys- 
tems are destroyed in the ensuing common envelope phase. 

The steeper Scalo IMF (Model F), normalised to the 
local space density of low-mass stars, produces fewer inter- 
mediate (and high) mass stars than the KTG IMF, and so 
fewer of the dominant progenitors in Fig. 14 are produced. 
More of the compact binaries are then descended from lower- 



Table 4. Summary of the properties of the optimistic and pes- 
simistic models, along with the fiducial Model A 



Model 


% DD 


(q) 


f2 gw (l mHz) 


(M) 


N 


Optimistic 


26 


0.75 


5.99 


0.46 


1.85 


A 


18 


0.62 


3.57 


0.45 


1.17 


Pessimistic 


14 


0.66 


0.95 


0.40 


0.32 



mass progenitors than for Model A, giving rise to their lower 
mean chirp mass. If we had instead normalised to the local 
core-collapse supernova rate, as in Schneider et al. (2001), we 
would instead have ended up with a correspondingly higher 
background from Model F. 

Altering the shape of the cosmic star formation history 
(Model J) has little impact upon the background, since most 
of the sources contributing have MS evolution times of less 
than a few Gyr (see Fig. 14). This is a strong argument in 
favour of using an integral constraint (such as IR luminosity 
density), and not a present-day constraint (such as local 
core-collapse supernova rate), since normalising according to 
the supernova rate introduces a strong dependence on the 
shape of the cosmic star formation history curve, through 
the difference in amplitude between the local rate and the 
rate at the peak of star formation, which can easily skew the 
overall normalisation. 

Finally, Models Q and R lead to larger gravitational 
wave backgrounds than Model A, mainly because lower 
metallicity stars tend to leave the main sequence earlier, and 
thus a greater fraction of the stellar mass in the universe to- 
day is present in the form of remnants. The difference in 
received flux is, however, slight, on the order of 10 per cent. 
We conclude that keeping detailed track of abundance vari- 
ations is not essential to calculation at the present level of 
accuracy. 



8 OUTLOOK 

Based on the above indications of which effects boost the 
GW background and which reduces it, we construct two 
models in an attempt to put upper and lower limits on the 
background we predict. Our use of the terms 'optimistic' and 
'pessimistic' assumes that this background constitutes signal 
for the reader; if it constitutes a noise, the nomenclature 
should be reversed. 

Optimistic model: This has the properties of Model 
A, except for: the Nelemans et al. (2000) common envelope 
formalism, initially circular orbits, enhanced mass transfer 
on the HG, edge lit detonations only after accretion of 0.3 
Mq and no ionisation energy in envelope binding energies 
used for common envelope phases. Note that some of these 
individually boosting effects do not make a double-boost in 
combination; for example the no spiral-in common envelope 
prescription tends to lead to DD mass ratios closer to unity, 
which means that fewer systems undergo stable mass trans- 
fer upon contact, and so the enhancement brought by the 
higher ELD limit is less effective in increasing the amplitude 
of the background. We also include the estimated error on 
our overall normalisation (see section 6.4), by using a cosmic 
star formation rate everywhere 30 per cent higher than our 
fiducial one. 
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Figure 16. Optimistic (upper dotted), fiducial (Model 
A, lower solid line) and pessimistic (lower dotted) extra- 
galactic backgrounds plotted against the LISA (dashed) 
single-arm Michclson combination sensitivity curve (see 
http://www.srl.caltech.edu/~shanc/sensitivity/). The 'unre- 
solved' Galactic close WD— WD spectrum from Nelcmans ct al. 
(2001c) is plotted (with signals from binaries resolved by LISA 
removed), as well as an extrapolated total, in which resolved 
binaries are restored, as well as an approximation to the Galactic 
MS-MS signal at low frequencies. 



Pessimistic model: The pessimistic model contains 
the elements found in the previous section to decrease the 
amplitude of the GW background. The properties of this 
model are thus the same as Model A, except for: a = 2.0, 
Roche lobe overflow is semiconservative on the HG, the Scalo 
initial mass function is used and 100 per cent of the ionisa- 
tion energy is included in the envelope binding energies used 
in common envelope phases. In addition, we use a star for- 
mation rate everywhere 30 per cent lower than our fiducial 
one, in our cosmological integral. 

These prescriptions were used to create Galactic DD 
populations, which were found to compare reasonably with 
observations. Then the cosmological integrals were carried 
out for each. The results of this are summarised in Table 4, 
and the optimistic, fiducial and pessimistic total background 
spectra are plotted in Fig. 16 along with the LISA sensitiv- 
ity curve, and the Galactic WD-WD background taken from 
Nelemans et al. (2001c). We plot both the 'unresolved' ('av- 
erage') background curve from their paper, which is for DD 
pairs only, with the resolved sources removed, and an ex- 
trapolated 'total' background. In this we have added back 
in the resolved close binaries and made an approximation to 
the MS-MS contribution at lower frequencies, in an attempt 
to represent the Galactic signal over the full frequency range 
plotted. 

Plotted in Fig. 17(a) is the number of systems per 1/(3 
yr) frequency resolution element contributing to the GW 
background as received today. We see from this that at fre- 
quencies / r < 50 mHz, there will be too many individual 
WD-WD sources contributing in each resolution element for 
this background to be completely resolved and subtracted 
source by source by missions with plausible lifetimes. How- 
ever, from Fig. 17(b), we see that much of the flux comes 
from relatively nearby sources, and the WD-WD numbers 
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Figure 17. (a) The number of systems per 10 Hz contributing 
to the cosmological GW background as received today. Linestyles 
denote the evolutionary classes as in Fig. 4. (b) Thin line: the 
fractional contribution at 10 mHz to the GW background as a 
function of cosmic time (from shells of width To/50). Thick line: 
the same, but in terms of the number of sources contributing to 
the flux received from each cosmological time-shell. 



drop rapidly above 50 mHz (leaving the lower background 
from rare neutron stars and black holes, not considered in 
this paper). Thus it may be possible for future missions more 
sensitive than LISA to subtract this background at high fre- 
quencies. 



9 COMPARISON WITH PREVIOUS WORK 

Hils et al. (1990) and Kosenko & Postnov (1998) each made 
an order of magnitude estimate of the ratio of the extra- 
galactic to Galactic GW flux from DDs. In order to facilitate 
comparison, and to compare like with like as far as possible, 
we divide our calculated extragalactic flux at 1 mHz by the 
most recently calculated value (Nelemans et al. 2001c) for 
the Galactic flux at the same frequency (which is a factor 
~ 3 smaller than that found by Webbink & Han 1998). The 
correct curve from Fig. 16 to use for this comparison is our 
'extrapolated' curve. We find f2 e xtra g ai /£l ga .i = 2.0 per cent 
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at 1 mHz for Model A, with a range of 0.5 - 3.4 per cent 
between optimistic and pessimistic Models. 

Hils et al. (1990) predicted a factor ~ 1.6 per cent (for 
an Einstein-de Sitter universe with no cosmological evolu- 
tion of galactic GW luminosity). This estimate is in good 
agreement with our value. 

Kosenko & Postnov (1998), on the other hand, pre- 
dicted that, for a cosmology of the type used in this pa- 
per, the extragalactic background should be of order 10 per 
cent of the Galactic one, when one takes into account the 
evolution of star formation rate with redshift. This result 
is in clear disagreement with our findings, but this can be 
explained by noting that their ratio is artificially raised by 
a number of factors: first, the fiducial scalings of fib, (r) 
and feioo in their eq. 13 are higher than their true values, 
boosting the extragalactic signal. Second, the same star for- 
mation rate as a function of redshift was used for different 
cosmologies, which leads to an artificial boost to the lambda- 
cosmology extragalactic flux (see e.g. Somerville, Primack & 
Faber 2001). Lastly, the cosmic star formation rate adopted 
was not normalised to any integral constraint, but merely to 
the current star formation rate. All of these factors lead to 
their calculation yielding a misleadingly high extragalactic 
contribution to the GW background. 

Schneider et al. (2001) made a more direct calculation 
of the background. At 1 mHz, their derived background level 
(for hioo — 0.7) is fi gw = 1.2 x 10~ n , with no quoted un- 
certainty on this value. This lies a factor two outside of our 
predicted range for the background. The discrepancy can be 
understood mainly in terms of their different method of nor- 
malisation: they normalised to the local core collapse super- 
nova rate, and used the steep Scalo IMF, meaning that more 
low- and intermediate-mass stars were born in their simu- 
lations than measured by Cole et al. (2001). As explained 
in Section 6.4, we believe that normalising to an integral 
constraint on the birth of low-mass stars is a more robust 
method. Schneider et al. (2001) also used a binary fraction 
of 100 per cent, cf. our 50 per cent. 

The shape of the spectrum in Schneider et al. (2001), 
however, we cannot explain. The spiral-in part of the spec- 
trum (/ r > 10~ 4 Hz) has the form expected from Section 4, 
but the static regime instead displays a prominent 'bump' 
at frequencies (/ r ~ 3 x 10 -5 Hz) just below the transi- 
tion to the spiral-in regime, the amplitude of which decays 
rapidly towards lower frequencies. No such feature is seen 
in our calculated spectra. This type of feature is difficult to 
explain in terms of the arguments in Section 4, unless the 
vast majority of WD-WD pairs are born precisely into this 
'bump', which seems unlikely, since the same feature is seen 
for all types of compact object pair (e.g. NS-NS, NS-BH), 
despite their very different formation routes. 

10 CONCLUSIONS 

We predict that the background of gravitational waves from 
extragalactic binary stars is 

(i) Dominated by double main sequence binaries for f T < 
10~ 4 Hz. 

(ii) Dominated by double white dwarf binaries for 10~ 4 < 

fr < 10 _1 HZ. 



Concentrating on the spectrum around 1 mHz: 

(i) The fraction of critical density in gravitational waves 
received in the logarithmic frequency interval around 1 mHz 
lies in the range 1 x 10~ 12 < fi gw < 6 x 10~ 12 , with the most 
likely value in the range 3 — 4 x 10~ 12 . 

(ii) The flux- weighted mean chirp mass of the contribut- 
ing binaries is (M) = 0.45 M©. 

(iii) Half of the background comes from binaries whose 
more massive (primary) star had a mass in the range 2-4 
M© (and ~70 per cent from primaries originally less mas- 
sive than 4 M©). The estimate of the background is thus 
more robust to uncertainties in the IMF and mass cuts if 
normalised to the present density of starlight than if nor- 
malised to core-collapse supernova rates. 

(iv) ~ 60 per cent of the GW signal is from binaries with 
initial semi-major axes in the range of 30-1000 stellar di- 
ameters, in which the Roche contact of both primary and 
secondary stars led to unstable transfer and a common en- 
velope. The background level produced by these systems 
is quite stable against uncertainties in the efficiency of the 
common envelope phase, though the signal can be changed 
somewhat through use of a non-standard common envelope 
prescription. 

(v) ~ 40 per cent of the GW flux comes from systems de- 
scended from binaries with initial semi-major axes of about 
5 stellar diameters, in which the first Roche contact occurred 
in the Hertzsprung gap, with stable overflow, but the sec- 
ond Roche contact led to unstable transfer and a common 
envelope. The background level produced by these systems 
is sensitive to uncertainties in common envelope and mass 
transfer physics. 

(vi) interacting systems (AM CVn binaries) contribute 
only about 10 per cent of the energy density in gravitational 
waves. 

The above holds true for 0.5 < / r (mHz) < 5. Above this 
range, as the lower-mass WD-WD pairs reach contact and 
drop out of the spectrum due to mergers, the properties 
change (values at 10 mHz in the parentheses which follow): 
the contribution from interacting binaries increases (26 per 
cent), the RLOF+CEE route contribution (44 per cent) and 
the mean primary progenitor mass increase (4.7 M©) and 
the mean chirp mass is higher (0.56 M©). 

We find that at all frequencies, our derived spectral 
shape can be understood in terms of simple arguments, and 
that this shape is essentially independent of the population 
synthesis model used. 
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